/*! 
 ***************************************************************************
 * \file BuckleyLeverett2D.cpp
 * Two phase flow in porous media in 2D.
 *
 * \par Description
 * This code solves the two phase flow in porous media in 1D, the
 * equations are posed in pressure-saturation formulation:
 * 
 * These equation are solved using the IMPES algorithm.
 * The domain of study is \f$ (x,y) \in [0,300 m] \times [0,1m] \f$, and is
 * initially filled with oil. Then, water is injected on the left side
 * and we study the traveling front between oil and water.
 *
 *  \par Compiling and running
 *  Modify the variables BASE and BLITZ in the file \c tuna-cfd-rules.in 
 *  according to your installation and then type the next commands:
 *  \verbatim 
   % make
   % ./BuckleyLeverett2D  \endverbatim
 ***************************************************************************
 *  \author    Luis M. de la Cruz [ Apr 17 16:17:33 GMT 2011 ]
 ***************************************************************************
 */

#include <iostream>
#include <cmath>
#include "Meshes/Uniform.hpp"
#include "Equations/TwoPhaseEquation.hpp"
#include "./BL/BLIP5.hpp"
#include "./BL/BLES5.hpp"
#include "Solvers/TDMA.hpp"
#include "Solvers/explicit.hpp"
#include "Utils/inout.hpp" 
#include "Utils/num_utils.hpp"
#include "Utils/GNUplot.hpp"
using namespace Tuna;
  
// Scalar field in two dimensions
typedef TunaArray<double,2>::huge ScalarField2D;

timer crono;

int main() 
{
  double length_x, length_y, dt, Tmax, tolerance = 1e-6; 
  int num_nodes_x, num_nodes_y, frequency, tdma_iter = 20;
  double porosity, permeability, mu_w, mu_o, Srw, Sro;
  double injection, pres_out, Sw_in;
  int idt_p;

  // Reading data from the input file
  std::ifstream input_cin ("inputBL2D"); 
  input_cin >> length_x >> length_y
	    >> num_nodes_x >> num_nodes_y
	    >> porosity
	    >> permeability
	    >> mu_w
	    >> mu_o
	    >> Srw
	    >> Sro
	    >> injection
	    >> pres_out
	    >> Sw_in
	    >> dt
	    >> Tmax
	    >> idt_p
	    >> frequency;
  input_cin.close();

    double Gamma = permeability;

    StructuredMesh<Uniform<double, 2> > mesh(length_x, num_nodes_x,
					     length_y, num_nodes_y);
    double dx = mesh.getDelta(X);
    double dy = mesh.getDelta(Y);
    mesh.print();

    ScalarField2D p ( mesh.getExtentVolumes() );
    ScalarField2D Sw ( mesh.getExtentVolumes() ); 
    ScalarField2D p_n( mesh.getExtentNodes()   );
    ScalarField2D Sw_n( mesh.getExtentNodes()   );

    // Initial conditions
    Sw = 0;
    p = pres_out; 
 
    Range all = Range::all();
    // Boundary condition
    Sw(0, all) = Sw_in;  

    NumUtils::interpolateToNodes(p_n, p);
    InOut::writeToFile_DX(p_n, 0, "./DataBL/2D/pres.", dx, dy);
    NumUtils::interpolateToNodes(Sw_n, Sw);
    InOut::writeToFile_DX(Sw_n, 0, "./DataBL/2D/satu.", dx, dy);
        
    SparseMatrix< Diagonal< double, 2> > A(num_nodes_x, num_nodes_y); 
    //    DiagonalMatrix< double, 2> A(num_nodes_x, num_nodes_y); 
    ScalarField2D              b(num_nodes_x, num_nodes_y);

    TwoPhaseEquation< BLIP5<double, 2> > pressure(p, A, b, mesh.getDeltas());
    pressure.setDeltaTime(dt); 
    pressure.setPermeability(permeability);
    pressure.setPorosity(porosity);
    pressure.setSrw(Srw);
    pressure.setSro(Sro);
    pressure.setViscosity_w(mu_w);
    pressure.setViscosity_o(mu_o);
    pressure.setNeumann  (LEFT_WALL, -injection * mu_w / permeability);
    pressure.setDirichlet(RIGHT_WALL);
    pressure.setNeumann  (TOP_WALL);
    pressure.setNeumann  (BOTTOM_WALL);
    pressure.setSaturation(Sw);
    pressure.print();

    //
    // In these case, the boundary conditions are calculated in
    // the numerical scheme (the adaptor BLES1).
    //
    TwoPhaseEquation< BLES5<double, 2> > saturation(Sw, A, b, mesh.getDeltas());
    saturation.setDeltaTime(dt); 
    saturation.setPermeability(permeability);
    saturation.setPorosity(porosity);
    saturation.setSrw(Srw);
    saturation.setSro(Sro);
    saturation.setViscosity_w(mu_w);  
    saturation.setViscosity_o(mu_o); 
    saturation.applyBounds(1, Sw.ubound(firstDim)-1);// just to move the indexes
    saturation.setNeumann(TOP_WALL);
    saturation.setNeumann(BOTTOM_WALL);
    saturation.setPressure(p);
    saturation.print();


    ofstream error_fp("./DataBL2D/error"), 
      residual_fp("./DataBL2D/residual");
    int iteration = 1, idprint = 1;
    double t = dt;
    double error = 1.0;
    double residual = 1.0;

    int pausa;

    int bi = saturation.get_bi(), bj = saturation.get_bj();
    int ei = saturation.get_ei(), ej = saturation.get_ej();

    frequency *= (86400 / dt);

    double total_time = 0, pres_time = 0, sat_time = 0;
    while (t <= Tmax) {

      if ( !(iteration % idt_p) ) {
	crono.tic();
	pressure.calcCoefficients();
	Solver::TDMA2DX(pressure, tolerance, tdma_iter, 1.0);
	pres_time += crono.toc();
	
	residual = pressure.calcResidual();
	error = pressure.calcErrorL1();
	
	error_fp << iteration << "\t" << error << std::endl;
	residual_fp << iteration << "\t" << residual << std::endl;	
	
	pressure.update();
      }  

      crono.tic();
      saturation.calcCoefficients();
      Solver::solExplicit2D(saturation);
      sat_time += crono.toc();

      error = saturation.calcErrorL1();
      saturation.update();


      /* */
      for(int j = 0; j <= (ej+1); ++j)
	for(int i = 0; i <= (ei+1); ++i) {
	  if ( Sw(i,j) > Sw_in) Sw(i,j) = Sw_in;
	  if ( Sw(i,j) < Srw) Sw(i,j) = Srw;	  
	}
      saturation.updatePhi(Sw);
      /* */
      
      if ( !(iteration % frequency) ) {
	NumUtils::interpolateToNodes(p_n, p);
	InOut::writeToFile_DX(p_n, idprint, "./DataBL2D/pres.",dx, dy);
	
	NumUtils::interpolateToNodes(Sw_n, Sw);
	InOut::writeToFile_DX(Sw_n, idprint, "./DataBL2D/satu.",dx, dy);
	std::cout << "\n\n ---- ID print => " << idprint
		  << "\t - Iteration => " << iteration
		  << "\t Step = " << t
		  << "\n ---- Res = " << residual
		  << "\t Err = " << error;
	std::cout.flush();
	++idprint;
      }

      ++iteration;
      t += dt;
    }

    cout << "\n\n Presure Elapsed time = " << pres_time << " segs \n";
    cout << "\n Saturation Elapsed time = " << sat_time << " segs \n";

    std::cout << "\n\n NOTE: The data were written in the [./DataBL2D] directory \n\n";

    return 0;
}




